close all
clear all
 
% Load RCP data and plot it
%load rf_rcp8_5.txt
load rf_rcp6_0_plus.txt
%load rf_rcp4_5.txt
%load rf_rcp2_6.txt
%yrs=rf_rcp8_5(:,1);
yrs=rf_rcp6_0_plus(:,1);
%yrs=rf_rcp4_5(:,1);
%yrs=rf_rcp2_6(:,1);
%forcing=rf_rcp8_5(:,3);
forcing=rf_rcp6_0_plus(:,3);
%forcing=rf_rcp4_5(:,3);
%forcing=rf_rcp2_6(:,3);
% Set start time in 2020
tstart=256;

%tend=336


%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify initial temperature anomaly in 2020
DeltaT0 =1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify model parameters %%%%%%%%%%%%%%%%%%
% Choose the mean reversion rate, Speed
% Speed is inverse of Hasselmann time
% lambda=1.2;   % Joules/metres^2/Kelvin
% Range of C quoted by Frame is 0.2 to 2 times 10^9 
C=1*1e9; % Watts/metre^2 /Kelvin 
% Choose noise variance 
% sigmaQ=0.25e8;
% sigmaQ=0
% Sigma=sigmaQ/C;

% We will measure time in years,
% Specify size of  time interval between points.
% so need a factor of number of seconds in year
% Effectively divides time in seconds by yearlength
dt=1; 
% So then have to multiply speed by yearlength
yearlength=3.1 *1e7;
%Speed = lambda*yearlength/C; % Kelvin per year 
% Display some useful quantities
%FirstYear=yrs(tstart)
%HasselmannTime=1/Speed
%phi=exp(-Speed*dt)

% And convert F/lambda to a function
%Level=ts2func(forcing(tstart:end)/lambda);

% Specify # of points at  which solution is wanted
nPeriods=length(forcing(tstart:end)); % Same as RCP dataset


% specify # of intermediate steps taken between the sample points
nSteps=12; % Monthly



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Do one evolution without noise 
Ntrials0=1;
% Sigma=0;
% Replace forming of the HWV object and call to HWV by a call to Hasselmann
%[DeltaTnoiseless,Times0]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials0,nPeriods,nSteps,dt);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Do many evolutions with noise
% Specify how many parallel worlds to simulate
% 800 works well
Ntrials=1000;
% Choose noise variance  
%  
sigmaQ=.9375e8;
% sigmaQ=0
%Sigma=sigmaQ/C;
% Replace forming of the HWV object   obj=hwv(Speed,Level,Sigma,'StartState',DeltaT0);
%and call to HWV [DeltaT,Times,Z]=simByEuler(obj,nPeriods,'Nsteps',nSteps,'Ntrials',Ntrials,'DeltaTime', dt);
% by a call to Hasselmann
%[DeltaT,Times]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials,nPeriods,nSteps,dt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ECS = 0.5:0.25:20
lambdas = 3.7 ./ ECS

parfor k = 1:1:length(lambdas)
    CS = ECS(k);
    Sigma=0;
    lambda=lambdas(k);   % Joules/metres^2/Kelvin 
    Speed = lambda*yearlength/C; % Kelvin per year 
    % And convert F/lambda to a function
    Level=ts2func(forcing(tstart:end)/lambda);

    [DeltaTnoiseless,Times0]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials0,nPeriods,nSteps,dt);

    Sigma=sigmaQ/C;
    [DeltaT,Times]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials,nPeriods,nSteps,dt);

    %Create a CSV write file
% First suppress the unit 3rd dimension in Delta T matrix
Temp=squeeze(DeltaT);
% Then concatenate 
M=[Times DeltaTnoiseless Temp];
 
filename = ['RCP60plusoutput_' num2str(CS) '.txt'];
csvwrite(filename,M)

end